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Abstract 

In [Ler53] and [Gar56], Leray and Carding have developed a multiplier technique for deriving 
a priori estimates for solutions to scalar hyperbolic equations in either the whole space or the torus. 
In particular, the arguments in [Ler53, Gar56] provide with at least one local multiplier and one local 
energy functional that is controlled along the evolution. The existence of such a local multiplier is 
the starting point of the argument by Rauch in [Rau72] for the derivation of semigroup estimates for 
hyperbolic initial boundary value problems. In this article, we explain how this multiplier technique can 
be adapted to the framework of finite difference approximations of transport equations. The technique 
applies to numerical schemes with arbitrarily many time levels, and encompasses a somehow magical 
trick that has been known for a long time for the leap-frog scheme. More importantly, the existence 
and properties of the local multiplier enable us to derive optimal semigroup estimates for fully discrete 
hyperbolic initial boundary value problems, which answers a problem raised by Trefethen, Kreiss 
and Wu [Tre84, KW93]. 


AMS classification: 65M06, 65M12, 35L03, 35L04. 
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Throughout this article, we use the notation 

'^:={CgC,|C|>1}, W:={CgC,\C\>1}, 

D ;= {C G C,|C| < 1}, ST={CgC,|C| = 1}, D:=DUS^ 

We let ^n(^) denote the set of n x A matrices with entries in K = M or C. If M G M* denotes 

the conjugate transpose of M. We let I denote the identity matrix or the identity operator when it acts 
on an infinite dimensional space. We use the same notation x* y for the Hermitian product of two vectors 
x,y € C” and for the Euclidean product of two vectors x,y G The norm of a vector x G C”’ is 
\x\ := (x* x)^/^. The induced matrix norm on is denoted || • ||. 

The letter C denotes a constant that may vary from line to line or within the same line. The dependence 
of the constants on the various parameters is made precise throughout the text. 
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In what follows, we let d > 1 denote a fixed integer, which will stand for the dimension of the space 
domain we are considering. We shall also use the space ^ of square integrable sequences. Sequences may 
be valued in for some integer k. Some sequences will be indexed by while some will be indexed by 
Z"^ or a subset of Z*^. We thus introduce some specihc notation for the norms. Let Ax* > 0 for i = 1,... ,d 
be d space steps. We shall make use of the £^(Z'^“^)-norm that we define as follows; for all v G £^(Z'^“^), 


h2(Zd-l) 


HaxA 




\k=2 / i=2 jiGZ 

The corresponding scalar product is denoted (•, Then for all integers mi < m 2 , we set 


m 2 




mi,m2 


;=Axi ^ 


\u 




Jl=mi 


to denote the .^^-norm on the set [mi, m 2 ] x Z'^"^ (mi may equal —00 and m 2 may equal + 00 ). The 
corresponding scalar product is denoted (•, Other notation throughout the text is meant to be 

self-explanatory. 


1 Introduction 

1.1 Some motivations and a brief reminder 

The ultimate goal of this article is to derive semigroup estimates for hnite difference approximations 
of hyperbolic initial boundary value problems. Up to now, the only available general stability theory 
for such numerical schemes is due to Gustafsson, Kreiss and Sundstrom [GKS72]. It relies on 
a Laplace transform with respect to the time variable, and the corresponding stability estimates are 
thereby restricted to zero initial data. A long standing problem in this line of research is, starting from 
the GKS stability estimates, which are resolvent type estimates, to incorporate nonzero initial data and 
to derive semigroup estimates, see, e.g., the discussion in [Tre84, section 4]. This problem is delicate for 
the following reason: the validity of the GKS stability estimate is known to be equivalent to a slightly 
stronger version of the resolvent estimate 

sup {\z\ - 1) ||( 2 ;I-r)“^||^(^ 2 (N)) < -koo, (I) 

where T is some bounded operator on ^^(N) that incorporates both the discretization of the hyperbolic 
equation and the numerical boundary conditions. Deriving an optimal semigroup estimate amounts to 
showing that T is power bounded. In finite dimension, the equivalence between power boundedness of T 
and the resolvent condition (1) is known as the Kreiss matrix Theorem, but the analogous equivalence 
is known to fail in general in infinite dimension. Worse, even the strong resolvent condition 

sup sup (jzj - 1)” WizI- < -koo , 

n>l 

does not imply in general that T is power bounded, see, e.g., the review [SW97] or [TE05] for details and 
historical comments. 

Optimal semigroup estimates have nevertheless been derived for some discretized hyperbolic initial 
boundary value problems. More specihcally, the first general derivation of semigroup estimates is due 
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to Wu [Wu95], whose analysis deals with numerical schemes with two time levels and scalar equations. 
The results in [Wu95] were extended by Gloria and the author in [CGll] to systems in arbitrary 
space dimension, but the arguments in [CGll] are still restricted to numerical schemes with two time 
levels. The present article gives, as far as we are aware of, the first systematic derivation of semigroup 
estimates for fully discrete hyperbolic initial boundary value problems in the case of numerical schemes 
with arbitrarily many time levels. It generalizes the arguments of [Wu95, CGll] and provides new insight 
for the construction of “dissipative” numerical boundary conditions for discretized evolution equations. 
Let us observe that the leap-frog scheme, with some specific boundary conditions, has been dealt with 
by Thomas [Tho72] by using a multiplier technique. It is precisely this technique which we aim at 
developing in a systematic fashion for numerical schemes with arbitrarily many time levels. In particular, 
we shall explain why the somehow magical multiplier + u'j for the leap-frog scheme, see, e.g., [RM67], 

follows from a general theory that is the analogue of the Leray-Carding method for partial differential 
equations, which we briefly recall now. 

The method by Leray and Carding [Ler53, Gar56] provides with suitable multipliers for scalar 
hyperbolic operators of arbitrary order. Namely, given an integer m > 0, we consider a partial differential 
operator of the form 

m+l 

k=l 

where t G M stands for the time variable, x G stands for the space variable^, and each operator Pk{dx) 
is a linear combination of spatial partial derivatives of order k: 

Pk{dx) = Pk,a dx > Qx '■= ■ ■ • dxt , \a\ := ai +■■■ Pad- 

\a\=k 


In the above formula, the Pk,aS are real numbers^. Well-posedness of the Cauchy problem 

Lu = 0, {u,dtu,...,d^u)\t=o = {uo,ui,...,Um), (2) 

in Sobolev spaces is known to be linked with hyperbolicity of L. Namely, if L is strictly hyperbolic, 
meaning that for all ^ G \ {0}, the (homogeneous) polynomial 

m+l 

Pk{iO:=i’^ ( 3 ) 

k=l \a\=k 

has m -|- I simple purely imaginary roots with respect to r, then the Cauchy problem (2) is well-posed in 
X • • • X In particular, there exists a constant C > 0, that is independent of the solution 

u and the initial data uq, ui,..., Um, such that there holds; 

m m 

sup ^ ||9f'u(t)||Hm-fe(Rd) < Cl X/ 111 LI™-'“(R'*) • (^) 

k=0 k=0 

The method by Leray and Carding gives a quick and elegant way to derive the estimate (4) assuming 
that the solution u to (2) is sufficiently smooth. By standard duality arguments, the validity of the a 

^The periodic case a; £ can be dealt with in a similar way and is actually the one considered in [Gar56]. 

^We restrict here for simplicity to linear operators with constant coefficients. 
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priori estimate (4) yields well-posedness -meaning existence, uniqueness and continuous dependence on 
the data- for (2). Hence the main point is to prove (4) assuming that u is sufficiently smooth and decaying 
at infinity so that all integration by parts arising in the computations are legitimate. The main idea is 
to find a suitable quantity M u, which we call a multiplier and that will be linear with respect to u, such 
that when integrating the quantity 0 = {Mu) {Lu) on the slab [0,r] x one gets the estimate (4) 
for free (negative times are obtained by changing t —?■ —t). Following [Ler53, Chapter VI] and [Gar56, 
Section 3], one possible choice of a multiplier is given by L' u where V stands for the partial differential 
operator of order m whose symbol is drP, with P given in (3). Why L' u is a good multiplier is justified in 
[Ler53, Gar56] . A well-known particular case is the choice of 2 dtu as a multiplier for the wave equation. 
Here P{t,^) = and therefore drP = 2r, hence the choice 2dtu. The latter quantity is indeed a 

suitable multiplier for the wave operator because of the formula^: 

d 

2dtu{dlu - A^u) = dt{{dtu)‘^ + Y^(dxjU)'^'j - 2 diVa; (dtu Va:u) . 


The important fact here is that the energy: 


d 

{dtu)‘^ + ^{dxjU)'^, 

is a positive definite quadratic form of the first order partial derivatives of u. Let us observe that the 
multiplier V u is local, meaning that its pointwise value at {t, x) only depends on u in a neighborhood of 
{t,x). This is important in view of using this multiplier in the study of initial boundary value problems. 
Another important remark is that the above energy is also local, and the arguments in [Ler53, Gar56] show 
that this property is not specific to the wave operator. The fact that both the multiplier and the energy 
are local is crucial in the arguments of [Rau72, Lemma 1]. In our framework of discretized equations, 
the multiplier will be local but the energy will not necessarily be so. We shall not exactly follow the 
arguments of [Rau72] which use time reversibility, but rather construct dissipative boundary conditions 
which will yield the optimal semigroup estimate we are aiming at. 

1.2 The main result 

We first set a few notations. We let Axi,..., Axd, At > 0 denote space and time steps where the 
ratios, the so-called Gourant-Friedrichs-Lewy parameters. A* := At/Axi, i = are fixed 

positive constants. We keep At £ (0,1] as a small parameter and let the space steps Axi,..., Axd vary 
accordingly. The £^-norms with respect to the space variables have been previously defined and thus 
depend on At and the GFL parameters through the mesh volume {Ax 2 • • • Axd on and Axi • • • Axd 

on Z'^). We always identify a sequence w indexed by either N (for time), Z'^“^ or Z'^ (for space), with 
the corresponding step function. In particular, we shall feel free to take Fourier or Laplace transforms of 
such sequences. 

For all j £ 7J^, we set j = {ji,j') with j' := {j 2 ,... ,jd) £ Z'^“^. We let p,q,r £ denote some fixed 
multi-integers, and define pi,qi,ri, p',q',r' according to the above notation. We also let s G N denote 

®We refer to [Gar56, page 74] for the generalization of such ’’integration by parts” formula to partial derivatives of higher 
order. 
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some fixed integer. We consider a recurrence relation of the form; 


s+l 

= At 

a=0 

s+1 


j G > 1, n > 0 . 


n+s+l 


+ = = l-ri,...,0, n>0, 


a=0 


.,n _ fn 

k ^3 ~ ^3 ’ 


where the operators and are given by: 


j G Z"', > 1 - n , n = 0,..., s . 


Pi 


<?1 9 




£i=—ri l'=—r' 


e^=0i'=-q' 


( 5 ) 


( 6 ) 


In (6), the a^^o-, ^e,ji ,<7 are rea/ numbers and are independent of the small parameter At (they may depend 
on the CFL parameters though), while S denotes the shift operator on the space grid: {S^v)j := Vj^£ for 
j, i G Z'^. We have also used the short notation 

p' d Pi q' d qi 

E -E E . E -E E ■ 

l'=-r' 1=2 li=-ri £'=-q' i=2 li=-qi 

The numerical scheme (5) is understood as follows: one starts with initial data (/j ), {fj) dehned 
for ji > 1 — ri- Assuming that the solution has been defined up to some time index n + s, n > 0, then the 
first and second equations in (5) should uniquely determine for ji > 1 — ri, j' G . The meshes 

associated with ji > 1 correspond to the interior domain while those associated with ji = 1 — ri,..., 0 
represent the discrete boundary. We wish to deal here simultaneously with explicit and implicit schemes 
and therefore make the following solvability assumption. 

Assumption 1 (Solvability of (5)). The operator Qs+i is an isomorphism on (.‘^{73). Moreover, for all 
{Fj) G f^(N* X Z‘^“^) and for all ..., go^. G £^(Z'^“^), there exists a unique solution (uj)j^>i_,.^ G 

to the equations 

f Qs+i Uj = Fj , j , ji > 1, 

\uj + Bj^^s+iui^j! = gj , j'GZ"*, ji = 1 - n, • • • ,0. 

In particular, Assumption 1 is trivially satisfied in the case of explicit schemes for which Qs+i is the 
identity (a£,s+i = ^ii,o ''' <^ 4,0 hr (6), with S the Kronecker symbol). 

The first and second equations in (5) therefore uniquely determine for ji > 1 — ri, and one 

then proceeds to the following time index n + s + 2. Existence and uniqueness of a solution {u'j) to (5) 
follows from Assumption 1, so the last requirement for well-posedness is continuous dependence of the 
solution on the three possible source terms {FJ^), {gj'), if^)- This is a stability problem for which several 
definitions can be chosen according to the functional framework. The following one dates back to [GKS72] 
in one space dimension and was also considered by Michelson [Mic83] in several space dimensions. It is 
specifically relevant when the boundary conditions are non-homogeneous (( 5 ”) ^ 0): 

Definition 1 (Strong stability). The finite difference approximation (5) is said to be ’’strongly stable” 
if there exists a constant C such that for all 'y > 0 and all At G (0,1], the solution (u’j) to (5) with 
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(/?) = • • • = [f?) = 0 satisfies the estimate: 


7 


7 At + 1 


E 


— 2'ynAt |||^n|||2 


lll-ri,+oo 


+ At( 


— 2'y n At 


n>sH-l 


n>s+l 


Pi 

E 

j=l-ri 


\U' 


n \\2 

ji, 


< c 


7 At + 1 
7 


^ 114,-lfe(»-i) 7 P) 


n>s+l 


n>s+l 


ji=l-ri 


The main contributions in [GKS72, Mic83] are to show that strong stability can be characterized by a 
certain algebraic condition, which is usually referred to as the Uniform Kreiss-Lopatinskii Condition, 
see [Coul3] for an overview of such results. We do not pursue such arguments here but rather assume 
from the start that (5) is strongly stable. We can thus control, with zero initial data, type norms of the 
solution to (5). Our goal is to understand which kind of stability estimate holds for the solution to (5) 
when one now considers nonzero initial data {fj), ■ ■ ■, (fj) in Our main assumption is the following. 

Assumption 2 (Stability for the discrete Cauchy problem). For all ^ G the dispersion relation 


s+l 

E 

(T = 0 


Qa[e 




Fid 


p 


= 0 , Qa{K) ■= Y ^ 


( 8 ) 


has s + 1 simple roots in D. (The von Neumann condition is said to hold when the roots are located in 
D.J In (8), we have used the elassical notation 


for K G (C\ {0})'^ and I G 

From Assumption 1, we know that Q<j+i is an isomorphism on , which implies by Fourier analysis 
that Qs+i(e®^i,...,e*^'*) does not vanish for any ^ G M'^. In particular, the dispersion relation (8) is a 
polynomial equation of degree s + 1 in z for any ^ G M'^. We now make the following assumption, which 
already appeared in [CKS72, Mic83] and several other works on the same topic. 

Assumption 3 (Noncharacteristic discrete boundary). For ii = —ri,... ,pi, z G C and rj G let us 

define 

s+i p' 

aeiiz,r]) :=Y^‘' . (9) 

£7=0 i'=—r' 

Then o-n and do not vanish on x and they have nonzero degree with respect to z for all 

T] G 

Our main result is comparable with [Wu95, Theorem 3.3] and [CCll, Theorems 2.4 and 3.5] and shows 
that strong stability (or ”CKS stability”) is a sufficient condition for incorporating 1“^ initial conditions 
in (5) and proving optimal semigroup estimates. The main price to pay in Assumption 2 is that the 
roots of the dispersion relation (8), which are nothing but the eigenvalues of the so-called amplification 
matrix for the Cauchy problem, need to be simple. This property is satisfied for instance by the leap-frog 
and modified leap-frog schemes in several space dimensions, under an appropriate CFL condition, see 
Paragraph 1.3. Our main result reads as follows. 


6 




Theorem 1. Let Assumptions 1, 2 and 3 he satisfied, and assume that the scheme (5) is strongly stable 
in the sense of Definition 1. Then there exists a constant C such that for all j > 0 and all At G (0,1], 
the solution to (5) satisfies the estimate: 


sup mu’"III? 


7 


n>0 


7Ai + l 




—2 7 n At IW^'n- 1|| 2 


II1—ri,+oo 


n>0 


+ ^ At( 


—2 7 n At 


Pi 


n>0 


ir“il,-172(2^-1) ^ 

jl=l-ri 


cr=0 


fcr|||2 I 7 1 ^ A. -27 n At III pn III 2 

111—ri,H-oo ' / j ^ lll-^ nil ,+00 


7 


n>s+l 

0 






• (10) 


n>s+l 




In particular, the scheme (5) is ’’semigroup stable” in the sense that there exists a constant C such that 
for all At G (0,1], the solution (n’t) to (5) with (F^) = (g^) = 0 satisfies the estimate 


S 


sup |||w”|||?_ 

n>0 


< 


CY.\inil-ru+oo- 


( 11 ) 


The scheme (5) is also -stable with respect to boundary data, see [Tre84, Definition f.S], in the sense 
that there exists a constant C such that for all At G (0,1], the solution (uj) to (5) with (F^) = (/”) = 0 
satisfies the estimate 

0 

sup < C Yj Yj Il5ii,-Il£2(z^-1)- 

Theorem 1 gives the optimal semigroup estimate (11), and is therefore an improvement with respect to 
our earlier work [Coul4] where in one space dimension, and under an appropriate non-glancing condition^, 
we were able to derive the estimate (here ri = r, pi = p since d = 1): 


7 


7 At + 1 




u' 


n 1 2 


n>0 


n>0 


j=l-r 


< 


E 


a|||2 , 7At + l 

III 1—r,+oo 

a=0 n>s+l 


0 


n>s+l j=l—r 

The latter estimate does not incorporate on the left hand side the quantity: 

sun p- 27 »iAt |||„,n|||2 

oup C III tX III _|_QQ 5 

n>0 


n\2 

j 


and was unfortunately still not sufficient for deriving the semigroup estimate (11). Our main contribution 
in this article is to exhibit a suitable multiplier for the multistep recurrence relation in (5). With this 
multiplier, we can readily show that, for zero initial data, the (discrete) derivative of an energy can be 
controlled, as in [Rau72], by the trace estimate of (up and this is where strong stability comes into play. 

^The non-glancing condition is unfortunately not met by the leap-frog scheme. 
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This first argument gives Theorem 1 for zero initial data (and even for nonzero initial data if the non¬ 
glancing condition of [Could] is satisfied). By linearity we can then reduce to the case of zero forcing terms 
in the interior and on the boundary. The next arguments in [Rau72] use time reversibility, which basically 
always fails for numerical schemes^. Hence we must find another argument for dealing with nonzero initial 
data. Hopefully, the properties of our multiplier enable us to construct an auxiliary problem, where we 
modify the boundary conditions of (5), and for which we can prove optimal semigroup and trace estimates 
by ’’hand-made” calculations. In other words, we exhibit an alternative set of boundary conditions that 
yields strict dissipativity. Using these auxiliary numerical boundary conditions, the proof of Theorem 
1 follows from a standard superposition argument, see, e.g., [BGS07, Section 4.5] for partial differential 
equations or [Wu95, CGll] for numerical schemes. 

Remark 1. Assumption 3 excludes the case of explicit two level schemes for which s = 0 and Qi = I, 
for in that case a-r^ and/or do not depend on z. However, this case has already been dealt with in 
[Wu95, CGll], and we shall see in Section 3 where the assumption that O-r^ and are not constant is 
involved, and why the proof is actually simpler in the case s = 0 and Qi = I. 


1.3 Examples 

1.3.1 One space dimension 

Our goal is to approximate the outgoing transport equation (d = 1 here): 

dtu + adxU = 0 , u\t=o = uo, (12) 

with t,x > 0 and a < 0. The latter transport equation does not require any boundary condition at 
X = 0. However, discretizing (12) usually requires prescribing numerical boundary conditions, unless one 
considers an upwind type scheme with a space stencil ”on the right” (meaning ri = 0 in (5)). We now 
detail two possible multistep schemes for discretizing (12). Both are obtained by the so-called method 
of lines, which amounts to first discretizing the space derivative dxU and then choosing an integration 
technique for discretizing the time evolution, see [GK095]. 


The leap-frog scheme. It is obtained by approximating the space derivative dxU by the centered 
difference {uj+i — Uj-i)/{2 Ax), and by then applying the so-called Nystrom method of order 2, see 
[HNW93, Ghapter HI.l]. The resulting approximation reads 

u]+^ + Xa{u]+l-u]/l)-u] = 0, 

which corresponds to s = p = r = 1. Recall that A > 0 denotes the fixed ratio At/Ax. Even though (12) 
does not require any boundary condition at x = 0, the leap-frog scheme stencil includes one point to the 
left, and we therefore need to prescribe some numerical boundary condition at j = 0. One possibility® 
is to prescribe the homogeneous or inhomogeneous Dirichlet boundary condition. With general source 
terms, the corresponding scheme reads 


' ^ ^ ^ (^n+l _ ^n+1^ _ ^ pn+2 ^ ^ ^ n>Q, 

n > 0 , (13) 

,//), J>0. 

®With the notable exception of the leap-frog scheme that is indeed time reversible ! 

®This is of course not the only possibility and we refer to [GK095, 01i74, Slo83, Tre84] for some other possible choices 
which might be more meaningful from a consistency and accuracy point of view. Our main concern here is a discussion on 
stability for (5) and the Dirichlet boundary conditions are a good illustration for this aspect. 


n+2 


— % ■ 


[{uU]) = {n 




Assumption 1 is trivially satisfied because (13) is explicit. The leap-frog scheme satisfies Assumption 2 
provided that A |o| < 1. In that case, the two roots to the dispersion relation 

2 ^ + 2 z A a sin ^2 — 1 = 0, 

are simple and have modulus 1 for all ^ G M. Assumption 3 is satisfied as long as the velocity a is nonzero, 
for in that case 01 ( 2 ;) = —a-i{z) = Xaz. The scheme (13) is known to be strongly stable, see [GT81]. 
In particular. Theorem 1 shows that (13) is semigroup stable. An illustration of this stability property 
is given in the numerical simulation of a bump function, propagating at speed o = — 1 towards the left. 
Homogeneous Dirichlet boundary conditions are enforced at j = 0. The reflection of the bump generates 
a highly oscillatory wave packet that propagates with velocity +1 towards the right. The envelope of this 
wave packet coincides with the prohle of the initial condition, which indicates that the £^-norm is roughly 
preserved by the evolution. This numerical observation is in agreement with semigroup boundedness. 

Other choices of numerical boundary conditions for the leap-frog scheme or its fourth order extension 
are discussed, e.g., in [01i74, Slo83, Tho72, Tre84]. The main discussion in [01i74, Slo83, Tre84] is to verify 
strong stability for a wide choice of numerical boundary conditions, and if strong stability holds, then 
Theorem 1 automatically gives semigroup boundedness, which was not achieved in these earlier works. 






Figure 1: Reflection of a bump by the leap-frog scheme with homogeneous Dirichlet condition at four 
successive times. 
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A scheme based on the backwards differentiation rule. We still start from the transport equation 
(12), approximate the space derivative dxU by the centered finite difference (ttj+i — Uj_i)/(2 Ax), and 
then apply the backwards differentiation formula of order 2, see [HNW93, Chapter III.l]. The resulting 
scheme reads; 


- H- 

2 ^ 2 




-2u 


n+l 




= 0 . 


This corresponds to s = 1 and 


Q 2 


(Uj)j6Z ^ Q Uj + ^ {Uj+l - Uj_l)^ 




The operator Q 2 is an isomorphism on f^(Z) since Q 2 is an isomorphism for any small A a (as a perturbation 
of 3/21), Q 2 depends continuously on A a, and there holds (uniformly with respect to A a): 


2 lll^lll—00,+00 — III Q 2 'u||| —oo, + oo • 

The operator Q 2 is therefore an isomorphism on £^(Z) for any Aa > 0 (see, e.g., [Cou09, Lemma 4.3]). 
Let us now study the dispersion relation ( 8 ), which reads here 


-+iAasin^j z^ —2z + - = 0. 

^ j Zi 


It is clear that the latter equation has two simple roots in z for any G M. Moreover, if sin^ = 0, the roots 
are 1 and 1/3 which belong to D. In the case sin^ 7 ^ 0, none of the roots belongs to and examining 
the case Ao sin^ = 1 , we find that for sin^ 7 ^ 0 , both roots belong to D (which is consistent with the 
shape of the stability region for the backwards differentiation formula of order 2, see [HW96, Chapter 
V.l]). Assumption 2 is therefore satisfied. Assumption 3 is satisfied as long as a is nonzero since there 
holds p = r = 1 and ai(z) = a_i(z) = Aaz^/ 2 . 

Theorem 1 therefore yields semigroup boundedness as long as one uses numerical boundary conditions 
for which the numerical scheme is well-defined (this is at least the case for A a small enough) and strong 
stability holds. 


1.3.2 Two space dimensions 


Here we wish to approximate the two-dimensional transport equation {d = 2): 

dtu + ai dx^u + a 2 dx^u = 0 , u\t=o = uq , 


in the space domain {xi > 0, X 2 G M}. When ai is negative, the latter problem does not necessitate any 
boundary condition at xi = 0. Following [AG76], we use one of the following two-dimensional versions of 
the leap-frog scheme, either 


.n+2 

^j,k 




n+l ' 

TA- 1 - 


— u 




= 0 , 


(14) 


or 


..n+2 , \ „ / ^j+l,fe+l + ^i+l,fc-l 

^j,k + “1 -- 2 --^ 


_L ,,"+1 


-|- A 2 02 


',,"■+1 - 77 ;"'+^ 

'^j+l,k+l ^ ‘^j-l,k+l 




-ulk = 0. 


(15) 
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Assumption 1 is trivially satisfied because (14) and (15) are explicit schemes. The scheme (14) satisfies 
Assumption 2 if and only if Ai |ai| + A 2 |a 2 | < I 5 while the scheme (15) satisfies Assumption 2 if and only 
if max(Ai |ai|, A 2 |a 2 |) < 1- Let us now study when Assumption 3 is valid. For the scheme (14), we have 
= Pi = 1, and 

ai{z, T]) = Xiaiz, a-i{z, rj) = -ai{z, rf ), 

so Assumption 3 is valid as long as oi 7 ^ 0. For the scheme (15), we have again ri = pi = 1, and 

ai{z, p) = z (Ai ai cosp + i A 2 02 sinr/), a_i(z, p) = z (—Ai oi cos 77 + i A 2 02 sinp ), 

so Assumption 3 is valid as long as both ai and 02 are nonzero. We refer to [AG79] for the verification 
of strong stability depending on the choice of some numerical boundary conditions for (14) or (15). Once 
again, if strong stability holds, then Theorem 1 yields semigroup boundedness and ^^-stability with respect 
to boundary data. 


2 The Leray-Garding method for fully discrete Cauchy problems 


This section is devoted to proving stability estimates for discretized Cauchy problems, which is the first 
step before considering the discretized initial boundary value problem (5). More precisely, we consider 
the simpler case of the whole space j G Z'^, and the recurrence relation: 


s+l 




i G Z^ 

n > 0 

cr =0 



II 


n = 0 


(16) 

.,s. 


where the operators Qa are given by (6). We recall that in (6), the are real numbers and are 
independent of the small parameter At (they may depend on the CFL parameters Ai,..., A^), while S 
denotes the shift operator on the space grid: {S^v)j := Vj+e for j,l G Z*^. Stability of (16) is defined as 
follows. 


Definition 2 (Stability for the discrete Cauchy problem). The numerical scheme defined by (16) is 
stable if Qs+i is an isomorphism from onto itself, and if furthermore there exists a constant Cq > 0 

such that for all At G (0,1], for all initial conditions {f^)j^zd ,..., {fj)j^id in there holds 


s 


sup 

nSN 


III 2 

III —00,+00 


a=0 


(17) 


Let us quickly recall that stability in the sense of Definition 2 is in fact independent of At G (0,1] (because 
(16) does not involve At and (17) can be simplified on either side by f[j Ax*), and can be characterized 
in terms of the uniform power boundedness of the so-called amplification matrix 


si/ (k) 


^ (^)/f5s-t-l (^) . 

1 0 .. 

0 

Vo 01 


-Qo(k)/Q6-hi(k)\ 

0 


G .^s+i(C), 


0 / 


(18) 


where the Qo-('^)’s are defined in (8) and where it is understood that is defined on the largest open 
set of C'^ on which Qs+i does not vanish. Let us also recall that if Qs-i-i is an isomorphism from 
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onto itself, then Qs+i does not vanish on and therefore does not vanish on an open neighborhood 

of With the above definition (18) for £/, the following well-known result holds: 

Proposition 1 (Characterization of stability for the fully discrete Cauchy problem). Assume that Qs+i 
is an isomorphism from ^^(Z'^) onto itself. Then the seheme (16) is stable in the sense of Definition 2 if 
and only if there exists a constant Ci > 0 such that the amplifieation matrix in (18) satisfies 


VnEN, VCeM'^, 





<Ci 


(19) 


In particular, the speetral radius of should not be larger than 1 (the so-ealled von Neu¬ 

mann condition). 

The eigenvalues of i 2 /(e*^i,..., e*^”^) are the roots to the dispersion relation (8). When these roots are 
simple for all ^ E the von Neumann condition is both necessary and sufficient for stability of (16), see, 
e. g., [Coul3, Proposition 3]. Assumption 2 is therefore a way to assume that (16) is stable for the discrete 
Cauchy problem. Our goal is to derive the semigroup estimate (17) not by applying Fourier transform to 
(16) and using uniform power boundedness of £/, but rather by multiplying the first equation in (16) by 
a suitable local multiplier. The analysis relies first on the simpler case where one only considers the time 
evolution and no additional space variable. 


2.1 Stable recurrence relations 

In this Paragraph, we consider sequences (u”')neN with values in C. The index n should be thought of as 
the discrete time variable, and we therefore introduce the new notation T for the shift operator on the 
time grid: := for all m, n E N. We start with the following elementary but crucial Lemma, 

which is the analogue of [Gar56, Lemme 1.1]. 

Lemma 1 (The energy-dissipation balance law). Let P E C[X] be a polynomial of degree s -|- 1 whose 
roots are simple and located in D. Then there exists a positive definite Hermitian form qe on and a 

nonnegative Hermitian form on that both depend in a way on P, such that for any sequence 

(r;'^),igN with values in C, there holds 

Vn E N , 2 Re(T(P'(T)u") P(T) u") = (s + 1) |P(T) + (T - /) {q^iv^,..., + qd{v^,..., u^+"). 

In particular, for all sequence that satisfies the recurrence relation 

VnEN, P(T)u’^ = 0, 

the sequence {qe{v'^,... ,v'^~^^))n£N is nonincreasing. 

The fact that there exists a Hermitian norm on that is nonincreasing along solutions to the 

recurrence relation is not new. In fact, it is easily seen to be a consequence of the Kreiss matrix Theorem, 
see [SW97]. However, the important point here is that we can construct a multiplier that yields directly 
the ’’energy boundedness” (or decay). The fact that the coefficients of this multiplier are integer multiples 
of the coefficients of P will be crucial in the analysis of Section 3, see also Proposition 2 below. 

Proof. We borrow some ideas from [Gar56, Lemme 1.1] and introduce the interpolation polynomials: 

yk = l,...,s + l, PkiX) := a - Xj), 


12 





where rci,... ,Xs+i denote the roots of P, and o 7 ^ 0 its dominant coefficient. Since the roots of P are 
pairwise distinct, the P^’s form a basis of Cs[X] and they depend in a way on the coefficients of P. 
We have 

s+l 

p' = Y.Pi^- 

k=l 

We then consider a sequence (r;"')„,eN with values in C and compute 
2 Re(T{P'{T)v^) P(T) v^^-{s + 1) |P(T) v^\‘^ 

S + 1 

= ^ T(Pfc(T))r;- (T - Xk) Pfc(T) v^ + T (Pfc(T) v^) (T - 5^) Pk{T)v- 

k=l 

s+1 

- - i^)(Pfc(T)r;-) (T - Xfc)(Pfc(T) v^) 

k=l 

s+1 

= J](T-|xfc|2)|P,(T)t;"|2. 

k=l 

The conclusion follows by defining: 


0 

> 

.,t(;®)EC^+\ qe{w^,-- 

s+l 

.,w^) := j;|Pfc(T)u;0|2, 

(20) 


qd{w^, ■ ■ 

k=l 

s+l 

.,w^) := J^(l-|xfcp)|Pfc(T)u;0|2. 
k=l 

(21) 


The form Qe is positive definite because the P^’s form a basis of Cs[W]. The form qa is nonnegative because 
the roots of P are located in D. Both forms depend in a way on the coefficients of P because the 
roots of P are simple. □ 

Lemma 1 shows that the polynomial P' yields the good multiplier TP'(T)u"' for the recurrence 
relation P(T) v'^ = 0. Of course, P' is not the only possible choice, though it will be our favorite one in 
what follows. As in [Gar56, Lemme 1.1], any polynomial of the form' 

s+l 

Q • ^ ^ ^k ^k ? ^1; • • • ? ^s+1 P d ? 

k=l 

provides with an energy balance of the form 

2Re (T(g(T)u-)P(T)u’^) =(«! + • • • + «,+!) |P(T)u"|2 + (T-/)(geK,...,^"+^)) + (?dK,...,^"+^), 

with suitable Hermitian forms qe,qd that have the same properties as stated in Lemma 1. 

'^The sign condition here on the coefficients a*, is the analogue of the separation condition for the roots in [Ler53, Gar56]. 
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2.2 The energy-dissipation balance for finite difference schemes 

In this Paragraph, we consider the numerical scheme (16). We introduce the following notation: 

S+1 S + 1 

( 22 ) 

a=0 cr=0 

Thanks to Fourier analysis, Lemma 1 easily gives the following result: 

Proposition 2 (The energy-dissipation balance law). Let Assumptions 1 and 2 he satisfied. Then 
there exist a eontinuous coercive quadratie form Eq and a continuous nonnegative quadratic form Dq 
on 1'^(Z‘^; M)^+^ such that for all sequences {v"')n^n with values in l'^(Z'^;M) and for all n G N, there holds 

2 (M L u”)_oo,+oo = (5 + 1) lll^ ^"lll'-oo,+oo + (T - /) Eoiv^, ..., u"+*) + Do{v^, u"+^) . 

In particular, for all initial data f^, ..., /^ G £^(Z‘^;M), the solution to (16) satisfies 

svipEo{v^,...,v^+^)<Eo{f,...,n, 

neN 


and (16) is (£‘^-)stable. 

Proof. We use the same notation v'^ for the sequence the corresponding step function on 

whose value on the cell [ji Axi, {ji + 1) Axi) x • • • x Ax^, {jd + 1) Ax^) equals x". Then Plancherel 
Theorem gives 

2 (Mx",Lx")_oo,+oo -{s + 1) |||i^;"||P-oo,+oo 

= 2 Re(T(P'(T))?(0) Pc(T) ^^(O) -(« + !) l^’c(T) ^^(OP ^ , 

where x” denotes the Fourier transform of x”, and where we have let 

s+l 

Pd^) ■= , Cj ■= , 

cr=0 

and P'i-{z) denotes the derivative of Pc_ with respect to z. 

From Assumption 2, we know that for all f G M'^, Pq has degree s + l and has s + l simple roots in 
D. We can apply Lemma 1 and get 


2 (Mx^Lx")_oo.+oo -{s + l) 

= [ {T-I)q,^^{^{0,...,i++^{0) + qd,d^i0,---,d^d0)7^, 

jRd (2vr)'^ 

where qe^(, qd,c_ depend in a way on () G and are 2 7r-periodic in each Q. Furthermore, q^^c, is positive 
definite and qd,c^ is nonnegative. The conclusion of Proposition 2 follows by a standard compactness 
argument for showing coercivity of Flo- D 
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2.3 Examples 

The first basic example corresponds to the case s = 0 for which the multiplier provided by Proposition 2 
is In that case, the energy Eo reads |||Qi'y|||?-oo,+oo (recall that Qi is an isomorphism) and the 

energy-dissipation balance law is nothing but the trivial identity 

2 (<3l «"+■ + Oo o”>-oo.+oo = IIIOl + Qo 

The second line of this algebraic identity can be rewritten as 

IIIQl t>”+'ll|ioo,+oo - IIIQl l>"lt-oo.+o„ + IIIQl t>”f-o„,+oo - IIIQo , 

and £^-stability for the Cauchy problem amounts to assuming that the operator norm of Qo is not 
larger than 1. Hence the dissipation term IQi _,_gQ — |||(5o r’”|||?.oo +oo nonnegative. 

Let us now consider the leap-frog scheme in one space dimension, for which we have s = 1 and 

L = T2 + AaT(S-S-i)-/. 

The corresponding dispersion relation (8) reduces to 

+ 2i \ a sin^ z — 1 = 0 . 

For A |a| < 1, the latter equation has two simple roots xi(^),X2(0 of modulus 1. Following the previous 
analysis, see (20)-(21), the form ge,c is given by 

— xi(C) -b — X 2 (C) = 2 -I- 2 -|- 4 AaRe(f sinCrc^ , 

and qd,(; is zero. The associated forms in Proposition 2 are Do = 0 and (recall here d = 1): 

Eo(u°, = 2 ^ Ax {Vjf + 2^Ax{vj)^ + 2Xa^Ax (u°+i - Vj_^) v] . 

j£Z j£Z j£Z 

The latter energy functional Eq is coercive under the condition A|o| < 1, which is the necessary and 
sufficient condition of stability for the leap-frog scheme, and Eq is conserved for solutions to the leap-frog 
scheme. The conservation of Eq is usually proved by starting from the recurrence relation 

VjeZ, VuGN, - Uj = 0 , 

using the multiplier + rtj, and summing with respect to j. This is equivalent, for solutions to the 
leap-frog scheme, to what we propose here, since our multiplier reads 

Mu] = 2 + A a (u"+i - u]+l) = + u] + L u] . 

=0 

However, it will appear more clearly in Section 3 why our choice for M u] has a major advantage when 
considering initial boundary value problems. 

Let us observe here that the energy functional Eq is associated with a local energy density 

Eojiv^, v^) := 2 + 2 {v]f + 2Xa (u^+i - v]_i) v] . 

This is very specific to the leap-frog scheme. In general, the coefficients of the Hermitian forms qe,(,Qd,( 
are not trigonometric polynomials of C and therefore Eq, Dq do not necessarily admit local densities. This 
is one main difference with [Ler53, Gar56]. 


15 


3 Semigroup estimates for fully discrete initial boundary value prob¬ 
lems 

We now turn to the proof of Theorem 1 for which we shall use the results of Section 2 as a toolbox. By 
linearity of (5), it is sufficient to prove Theorem 1 separately in the case (/j*) = • • • = (/?) = 0, and in the 
case (F^) = 0, (g^) = 0. The latter case is the most difficult and requires the introduction of an auxiliary 
set of “dissipative” boundary conditions. Solutions to (5) are always assumed to be real valued, which 
means that the data are real valued. For complex valued initial data and/or forcing terms, one just uses 
the linearity of (5). 


3.1 The case with zero initial data 

We first assume (/°) = • • • = [fj) = 0. By strong stability, we already know that (7) holds with a constant 
C that is independent of 7 > 0 and At G (0,1]. Therefore, proving Theorem 1 amounts to showing the 
existence of a constant C, that is independent of 7 > 0 and At G (0,1] such that the solution to (5) with 
(/?) = • • • = if-) = 0 satisfies 


sup e 

n>0 


—2'ynAt |||^n|||2 


lli_.„+oo < C E Ate-2^"^‘|||FTI|2 


7 


n>s-\-l 


l,+oo 




n>s+l 


jl=l-ri 


(23) 


We thus consider a parameter 7 > 0 and a time step At G (0,1], and focus on the numerical scheme (5) 
with zero initial data (that is, (/?) = • • • = (/?) = 0). For all n G N, we extend the sequence (u”) by zero 
for ji < -rf. 

„ J if ji > 1 - ri, 

Vj := < 

0 otherwise. 


We use Proposition 2 and compute: 


(T - I) Eo{v^, ..., + Do{v^, ..., = 2 (M L u-)_oo,+oo - (^ + 1) III III loo, +00 • 

Due to the form of the operator L, see (22), and the fact that vanishes for ji < —ri, there holds: 




if ji < -n - Pi , 


and we thus get 


(T - I) Eo{v^, ..., u"+*) + Do(u",..., u"+^) 

= ("p Axk\ E E 2 At (M v^) - (s + 1) At^ (ir«+^+i)2 

Va:=1 / ji>l 

+ E E 2{MvJ)LvJ-{,+ 1){Lv^)\ 

\k=l / ii=l-r-i-pi j'gZd-l 
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We multiply the latter equality by exp(—2 7 (n + s + 1) At), sum with respect to n from 0 to some N and 
use the fact that Dq is nonnegative. Recalling that the initial data in (5) vanish, we get 


N 


^-2^{N+s+l)At _ _ _ ^^N+s+1-^ + (1 _e- 27 At^ (n+s) At ^ ^ 


n=l 


>0 


< 5'i,iv + S 2 ,n , (24) 


with 


and 


N 

.= (2At(Mu’*,F’^+*+^)i,+oo - (s + 1) At^ , (25) 


n=0 

d 


\ Af 0 

S 2 ,Ar:=fn^^fc E E ‘^iMv^)Lv^-{s + l){Lv^f. (26) 

\fc=l / n=0 ji=l—fi—Pi j'eZ'*”! 

Let us now estimate the two source terms Si^n, S 2 ,n in (24). We begin with the term S 2 ,n defined in 
(26). Let us recall that the ratio At/Axi is fixed. Furthermore, the form of the operators L and M in 
( 22 ) gives the estimate (recall that u” vanishes for ji < —ri): 


s2,N<cAt E® 


N Pi 

—2'y (n+s+1) At 


\k=2 


n=0 




U 


n+s+l\2 


for a constant C that does not depend on N, 7 nor on At. We thus have, uniformly with respect to 
A G N, 7 > 0 and At G (0,1]: 


W+s+l pi 

—2 7 n At 


S2,n<C Ate 

n=s+l 


E 


\u 


n \\2 




<C Y 


Pl 


— 2') n At 


< c 


7 At + 1 

7 


jl=l-ri 

^—2 7 n At I 


n>sH-l 


^ Ate-2^"^*|||F-|||2_+^+ At 


n>5+l 


E Il'“ii,-llr2(zd-i) 

jl=l-ri 

^-2,nA, ^ ||j,„ I 

n>s+l ji=^—Ti 


d-l\ ( 5 


where we have used the trace estimate (7) that follows from the strong stability assumption. 

Let us now focus on the term in (24), see the defining equation (25). We use the Cauchy- 
Schwarz inequality and derive (using now the interior estimate in (7) that follows from the strong stability 
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assumption): 


N 


Si,N < 2 ^ Ate- 2 |||Mu"|||i,+oo |||F”+*+^|| 


l,+oo 


n=0 

N 


<CY, Ate-2^("+^+i)^* (|||u”+i|||i_,,,+oo + • • • + 


n=0 


||,,n+s+l||| 

11 *^ ll|l—?’l,+oo 


7in+s+l|| 


II l,+oo 


< c 


N+s+1 

'y \ '• A j. 27 n At III n|||2 


7 At + 


- ^ Aie-^”“|||u”|||L„.+^ + c 


7 At + 1 


N+s+1 


< c 


n=s+l 

7 At + 1 


7 




+00 


n=s+l 


7 


^ Ate-2^-^*|||F"|||2_+^+ Ate-2^-^‘ ^ \\gl^.\ 


n>5+l 


n>s+l 


jl=l-ri 


Ignoring the nonnegative term on the left hand-side of (24) and using the coercivity of Eq, we have proved 
that there exists a constant C > 0 that is uniform with respect to iV, 7 , At such that: 


,-27(iV+s+l)At|||^7V+s+l|||^^^^^ < C I ^ At< 

y ^ n>s+l 


— 2'ynAt in rpn\\\2 


II l,+oo 


0 


+ ^ llaS,. 


— 2'ynAt \ ^ 11^?^ \\2 

n>s+l ii=l—?"i 


-1\ / 5 


which yields (23) and therefore the validity of Theorem 1 in the case of zero initial data. 


3.2 Construction of dissipative boundary conditions 

In this paragraph, we consider an auxiliary problem for which we shall be able to prove simultaneously 
an optimal semigroup estimate and a trace estimate for the solution. More precisely, we shall prove the 
following result. 


Theorem 2. Let Assumptions 1, 2 and 3 he satisfied. Then for all Pi G N, there exists a constant 
Cp.^ > 0 such that, for all initial data {fj), ■ ■ ■, (fj) G and for all source term (g(7)jj<o^„>s+i that 

satisfies 

Vr > 0 , ^ e- 2 r»^ ^ < +00, 

n>s+l 

there exists a unique sequence solution to 


'Lu^ = 0, J G Zi'^, 

<Mu] = gY^\ 
u] = fy, jGZ'', 


ji > 1, n > 0, 
ji < 0 , n > 0 , 
n = 0,..., s . 


(27) 
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Moreover for all 'j > 0 and At G (0,1], this solution satisfies 


sup e 

n>0 


— 2^71 At ||L,'^|||2 


II —oo,+oo 


+ 


7 


7 At + 


I 


—2'ynAt |||^n|||2 


Pi 




n>0 


—2jnAt 
n>0 ji=l—ri 


7/” 11^ 


< 


cp, \ E iiiri-o<.,+o<. + E 

iT=0 n>s+l 


E iisJ.i 

ii<o 




• ( 28 ) 


Theorem 2 justifies why we advocate the choice Mu” = + Aa rather than the 

more standard as a multiplier for the leap-frog scheme. Despite repeated efforts, we have not 

been able to prove the estimate (28) when using the numerical boundary condition + vd- on ji < 0 , 
in conjunction with the leap-frog scheme on > 1 . 


Proof. Let us first quickly observe that the solution to (27) is well-defined since, as long as we have 
determined the solution up to a time index n -|- s, n > 0 , then is sought as a solution to an 

equation of the form 

= F, 

where F belongs to (this is due to the form of L and M, see (22)). Hence is uniquely dehned 

and belongs to for all n G N. 

The proof of Theorem 2 starts again with the application of Proposition 2. Using the nonnegativity 
of the dissipation form Dq, we get® 


(T - /) ^;o(n^ ... ,u"+^) + (s + 1) < 2 (Mu^Ln’^)_oo,+oo = 2 Ln")_oo,o • 


By the Young inequality, we get 

(T - I) Eo{u", ... ,u”+*) + ll|L«”ll|i»o,+» < jY „ . 

We multiply the latter inequality by exp (—27 (n -|- s -|- 1) At), sum from n = 0 to some arbitrary N and 
already derive the estimate (here we use again the fact that At/Axi is a fixed positive constant): 


sup e- 2 T'(’^+^) u’"+") + (1 - e-^T-^*) ^ ^tEo{u^,..., u’"+") 

n>0 


+ '^At( 


-2j{n+s+l)At 


Eiii 


||2 


n>0 


iisz 




r2(z‘'-i) 


n>s-|-l 


ii<0 


Using the coercivity of Eq and the inequality 


1 — e 


-2'yAty 7 At 


7 At -|- 1 ’ 


®Smce Luf = 0 for ji > 1, one could also write |||I/u’*|||?.oo,o rather than |||Lu’^|||?.oo,+oo on the left hand-side of the 
inequality. 
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we have therefore derived the estimate 


sup e 

n>0 


—2771 At |||^n|||2 


+ 


7 


II- 00,+00 ' 

+ ^ At 


I 


^ — 2777 At III ^77 III 2 


II —cxd,H-oo 


77>0 
— 27 (ti+s+I) At 


u 


n ||2 


7 


n>0 


< 


He 


Eiii 

+ ^ ^t^-2'rnAt 


'"|||^oo,+oo + > Ate 

(T = 0 77>S+1 


E ^ , (29) 

ii<o 


where the constant C is independent of 7 , At and on the solution (u^)- In order to prove (28), the main 
remaining task is to derive the trace estimate for (u7). This is done by first dealing with the case where 
7 At is large. 


From the definition of the operator L, see (22), there exists a constant C > 0 and an integer J such 


that 


(i u]? > ((Q.+17+*+')" - c H E (“i 


72+(J\2 

’j+i) 


cr=0 |£|<J 

Since Qs+i is an isomorphism, there exists a constant c > 0 such that 

1 


E ’‘If 2 c 5: - T E E (“ 7 ") 


o-=o jgzd 

Multiplying by exp (—27 (n + s + 1) At) and summing with respect to n G N, we get 


Ate-^T'^^At 


b/” iP ^ 


1) < C <j ^Ate-2T'(’^+"+i)^* ^ ||L 

n>0 JiSZ 


n>s+l 


iiez 


||2 


^e-27Ai ^^^g-27nAt ^ | 
n>0 ilSZ 


7 /” ||2 


Choosing 7 At large enough, that is 7 At > In Rq for some numerical constant i?o > 1 that depends only 
on the (fixed) coefficients of the operator L, we have derived the estimate 


Y Ate-^T'"^* Y 


w 


n ||2 




1 ) < C i ^ Ate-^T'^ ||L 


u 


n ||2 


n>s+l 


iisz 


n>0 


il6Z 


^g -2 7 At ^ 


p(j III 2 

' III—00,+00 


(7 = 0 


It remains to use (29) and we get an even better estimate than (28) which we were originally aiming at: 


Y Ate-^T'^At Y 


I'lt?, .Il*2/'7d-l^ E C* < ^ ^ 




72>S+1 




(7 III 2 

III —00,+00 

(7=0 77.>S + 1 


ii<o 


20 



This gives a control of infinitely many traces and not only finitely many (this restriction to finitely many 
traces will appear in the regime where 7 At can be small). 

• From now on, we have fixed a constant > 1 such that (28) holds for 7 At > lni?o and we thus 
assume 'j At G (0,lni?o]- We also know that the estimate (29) holds, independently of the value of 7 At, 
and we now wish to estimate the traces of the solution (uj) for finitely many values of ji. 

We first observe from (29) that for all 7 > 0 and At G (0,1], there exists a constant such that 

Vn G N , • 

jeZ'* 

In particular, for any ji G Z, the Laplace-Fourier transforms ilj) of the step functions 

d 

Uj, : (t,y) G M+X u" if (t, y) G [n At, (n + 1) At) x Ax^, (jfc + 1) Ax^), 

k=2 

is well-defined on {r G C, Re r > 0} x The dual variables are denoted r = 7 + i0, 7>0, and 

V = {V 2 , ■ ■ ■, rjd) £ It will also be convenient to introduce the notation r/A := (72 ^3^2; ■ ■ ■ ,Vd ^Xd)- 

Given F > 0, the sequence (ft 7 (F -|- i6,r]))j-^^z belongs to £^(Z) for almost every (0,r/) G M x 
We first show the following estimate. 

Lemma 2. With Rq > 1 fixed as above, there exists a constant C > 0 such that for all 'j > 0 and 
At G (0,1] satisfying j At G (0,lnRo]; there holds 


E 

ii 6 Z- 


Pi 




£i=-ri 


d9 dr] 


+ E 

ii<o- 


pi 


Y ^a) (7 + ^ 61 , 7 ) 

il=-ri 


d6 di] 


< 


c 1 E iiirP-o„,+o„ + E 


fo -|||2 

II —00,+00 

(T=0 n>s+l 


— 2'y n At 


ii<o 


ti=-ri 


cr=l (j'=0 

_ , _-TAt ^+1 ^-1 

^ e"^*5,a^,(e"^*,7?A)5^(r,7) = Mn,-,,(r,7) +-^- E (32) 


Pi 


T 

1 — e“ 


(30) 


Proof of Lemma 2. Given t = 'y + i 6 and 7 , we compute (here ji G Z is fixed): 

Pi 1 _p-TAt «+i ^-1 , 

E ae,{e^^\r]A) Sj^fiT,r]) = Luj,,iT,v) + --- EE e(‘^-"')"^'^E(77), (31) 


£l=-ri 

where, in (31) and (32), we have set 


(T=l cr'=0 


Pi 


^r'('rt= E E 1 .(')). 

£i=—ri \i'=—r' 
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which corresponds to the partial Fourier transform with respect to y = (x 2 , ■ ■ ■ ,Xd) € of the step 

function associated with the sequence (Qo- ff ) (no Laplace transform here). 

We need to estimate integrals with respect to (0,^) of the right hand side of (31) and (32). The first 
term on the right of (31) and (32) are easy. For instance, we have (applying Plancherel Theorem): 

'"(n+l) At 


^ / jLuj,,.(T,?])fd0d7] = (2 7r)'^ ^ ^ 

1 jiez n>0 




= ( 2 vr)^ 


jieZ n>0'"‘ 
1 -e-^^^* 

2 7 At 


At 


^Ate-2^"^* ^ \\L 


u 


L' 


n ||2 
'Jl,- 


n >0 iiSZ 

We now recall that 7 At is restricted to the interval (0,lni?o]) and we use (29) to derive 
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Similarly, we have 
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which we can again uniformly estimate by the right hand side of (30). 

Going back to the right hand side terms in (31) and (32), we find that there only remains for proving 
(30) to estimate the integral (here there are finitely many values of a and a')\ 


E 

iisz’ 


1 — e 


— T At 


\-^jf'{ri)\^dedr] = 


1 — e 


— T At 




where we have applied Fubini Theorem. Applying first Plancherel Theorem with respect to the d — 1 last 
space variables, we get 
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The conclusion then follows by computing 

1 — e 
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Jr 


— T At 2 .. _ —27 At 

d9 = 27r AT 


2 7 At ’ 

and by recalling that 7 At belongs to (0,lni2o]- We can eventually bound the integrals on the left hand 
side of (30) by estimating separately the integrals of each term on the right hand side of (31) and (32). □ 

The conclusion now relies on the following crucial result. 

Lemma 3 (The trace estimate). Let Assumptions 1, 2 and 3 he satisfied. Let Rq > 1 be fixed as above 
and let Pi G N. Then there exists a eonstant Cp^ > 0 such that for all z G ^ with \z\ < Rq, for all 
T] G and for all sequence (utji)jiez £ ^^(^; C), there holds 


(33) 




pi 

2 

pi 


E KlTCft. 

jl=-ri-pi 


£l=-ri 

+ E 

ii<o 

^ zdza£fiz,r])wj^+i^ 

il=-ri 



Recall that the functions a^j, £1 = —n ,... ,pi, are defined in (9). 
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The proof of Lemma 3 is rather long. Before giving it in full details, we indicate how Lemma 3 
yields the result of Theorem 2. We apply Lemma 3 to z = exp(T At), r = 7 + i 0 with 7 At G (0, lni?o]) 
7 A £ and the sequence (t, r/))j^g 2 - We then integrate (33) with respect to {9,r]) and use Lemma 
2 to derive 

(T=0 n>s+l 

It remains to apply Plancherel Theorem and we get 
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Recalling that 7 At is restricted to the interval (0, InRo]) we have thus derived the trace estimate 


Ea*' 

nGN 


—2 7 n At 


Pi 

E 

ji=-ri-pi 


U 




d- 1 ) < C 


s 

El 

(T = 0 


II —00,+ 00 


+ 


n>s+l 
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Combined with the semigroup and interior estimate (29), this gives the estimate (28) of Theorem 2 for 
7 At G (0, In Rq]. □ 


Proof of Lemma 3. Let us recall that the functions are 2 vr-periodic with respect to each coordinate of 
r/. We can therefore restrict to r/ G [0, 2 7 r]'^“^ rather than considering 7 G We argue by contradiction 

and assume that the conclusion to Lemma 3 does not hold. This means the following, up to normalizing 
and extracting subsequences; there exist three sequences (indexed by /c G N): 

• a sequence with values in t'^(Z;C) such that ... ,u>p7 belongs to the unit sphere 

of C^i+^i+Pi+i for all k, and ... ,Wp^) converges towards as k tends to 

infinity, 

• a sequence {z^)keN with values in ^ n {C G C, |C| < Ro}, which converges towards z £ , 

• a sequence {r]^)k&N with values in [0,2 7r]'^“^, which converges towards rj G [0,2 7r]^“^, 
and these sequences satisfy: 


lim 

fc^+OO 



pi 

2 

pi 

E 


+ E 

Y ^''dzae,{z’^,ri’^)wYei 

iiez 

£i=-ri 

ii<o 

ii=-ri 


(34) 


We are going to show that (34) implies that ... ,Wp^) must be zero, which will yield a contra¬ 

diction since this vector must have norm 1 . 
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• Let us first show that each component ji G Z, has a limit as k tends to infinity. This is 

already clear for ji = —ri — pi,..., Pi. For ji > Pi, we argue by induction. From (34), we have 

Pi 

£l=-ri 

and by Assumption 3, we know that ap^{z,r]) is nonzero. Hence (iPpj+i)fceN converges towards 

^ pi-i 

^ E^^«^i(^’!Z)^pi-pi+i+l > 

which we define as We can argue by induction in the same way for all indices ji > Pi + 1, but 

also for indices ji < —ri — pi because the function a-r^ also does not vanish on x 

Using (34), we have thus shown that for each ji G Z, (tc^JfceN tends towards some limit Wj.^ as k 
tends to infinity, and the sequence w, which does not necessarily belong to ^^(Z; C), satisfies the induction 
relations: 


Pi 

VjiGZ, «Lfe!/)%i+ri = 0> (35) 

£l=-ri 

Pl 

Vji<0, = 0. (36) 

£l=-ri 

• The induction relation (35) is the one that arises in [GKS72, Mic83] and all the works that deal 
with strong stability. The main novelty here is to use simultaneously (35) for controlling the unstable 
components of ... ,tP-i) and (36) for controlling the stable components of ... ,W-i)- 

The fact that w satishes simultaneously (35) and (36) for ji < 0 automatically annihilates the central 
components. This sketch of proof is made precise below. 

We define the source terms: 

pi Pl 

£l=—ri £i=—ri 


which, according to (34), satisfy 

lim y|p((|2 = o, lim y|G^|2 = 0. (37) 

iisz ii<o 

We also introduce the vectors (here T denotes transposition) 
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and the matrices: 


L(z,??) : = 


(-api-i{z,r])/ap^{z,ri) . -a-rAz,'!) / ap^{z,'q)\ 

1 0 ... 0 


V 


0 

0 


0 1 


0 


G (^) 


(38) 


/ 


(z,??) : = 


( dzapi-Az^r])/dzapAz^v) . dzCL—n {z, t])/dzCLp^ {z, Tj)'^ 

1 0 ... 0 


0 

0 


0 1 


G ^p-^^+rl{C■) . (39) 


The matrix L is well-defined on ^ x according to Assumption 3. The matrix M is also well-defined 
on X because for any rj G Assumption 3 asserts that a,pi{-,'r]) is a nonconstant polynomial 

whose roots lie in D. From the Gauss-Lucas Theorem, the roots of dzO.pA'^'n) li® i^^ fh® convex hull of 
those of apA--,fi)- Therefore dzapA--,^) does not vanish on In the same way, dza-rA'iV) does not 
vanish on ^. 

With our above notation, the vectors , satisfy the one step induction relations: 

VjiGZ, Wj;+l=L(z^r?")W/^ + (T,';+l/ap,(z^^7'^),0,...,0)^, = L(z, ^) , (40) 

Vii < -1, =M(z^r,")IF,^ + (^Gl^Aiz'^dzapAzAq’^)),0,...,oY , = M(z, 5 ) . 

(41) 


• From Assumption 3 and the above application of the Gauss-Lucas Theorem, we already know that 
both matrices l^{z, rj) and M(z, rj) are invertible for {z, rj) £ ^ x Furthermore, Assumption 2 shows 

that 'L{z,r]) has no eigenvalue on for (z,r/) G x This property dates back at least to [Kre68]. 

However, central eigenvalues on may occur for L when z belongs to S^. The crucial point for proving 
Lemma 3 is that Assumption 2 precludes central eigenvalues of M for all z G '2^. Namely, for all z G 
and all rj G M(z,r 7 ) has no eigenvalue on S^. This property holds because otherwise, for some 

(z, 77 ) G X there would exist a solution ki G to the dispersion relation 

Pi 

zdzae^{z,r]) = 0 . 

il=-ri 

For convenience, the coordinates of rj are denoted ( 772 ,..., rjA). Using the definition (9) of , and defining 
K := (ki, e*'^^,...,e*'^'^), we have found a root z G ^ to the relation 

^f+i 

^cr(3^(K)z‘^"^ = 0, (42) 

(7=1 

but this is not possible because the s -|- 1 roots (in z) to the dispersion relation (8) are simple and belong 
to D. The Gauss-Lucas Theorem thus shows that the roots to the relation (42) belong to D (and therefore 
not to '^). 

At this stage, we know that the eigenvalues of M(z, 77 ), (z, 77 ) G x split into two groups: those 

in which we call the unstable ones, and those in D, which we call the stable ones. For (z, 77 ) G 
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we then introduce the spectral projector ry), resp. ry), of M( 2 ;, rj) on the generalized eigenspace 

associated with eigenvalues in D, resp. ^. We can then integrate the first induction relation in (41) and 
get 


^^(z^ry"X 


1 


ii<o 



The projector depends continuously on {z,r]) G x Furthermore, since the spectrum of M 

does not meet even for z G there exists a constant C > 0 and a 5 G (0,1) that are independent of 
G N and such that 


Vji < 0, ||M(z^ry^)l^ll^^(z^ry^)|| < . 


We thus get a uniform estimate with respect to k: 


UUz\r^^)W^\‘^<CY,\Gjf 

ji<0 


Passing to the limit and using (37), we get n^(z, ry) Wn = 0, or in other words Wn = n^(z, ry) Wn - 

• The sequence (E()ji)ii<o satisfies both induction relations (40) and (41). Due to the form of the 
companion matrices L and M, see (38)-(39), we can conclude that the vector Wn belongs to the generalized 
eigenspace (of either L or M) associated with the common eigenvalues of M(z, ?y) and L( 2 , ry). We have 
already seen that M( 2 , ry) has no eigenvalue on and Wq = n^(z, ry)lTo, so we can conclude that 
VF n belongs to the generalized eigenspace of L associated with those common eigenvalues of M(, 2 , ry) and 
L(z, ry) in 

The matrix L(z, ry) has N'^ eigenvalues in in D and iV^ on §^. (Since z may belong to 

N'^ is not necessarily zero.) With obvious notations, we let n“’*’'^(z, ry) denote the corresponding 
spectral projectors of L for (z,ry) sufficiently close to ( 2 , ry). In particular, the eigenvalues corresponding 
to nj(( 2 ;,ry) lie in uniformly away from for (^,ry) sufficiently close to ( 2 ,ry). We can then integrate 
the first induction relation in (40) and derive (for k sufficiently large): 


Iil{z\T,^)W^ = 


1 


j;L(/,ry^ 

ii>o 



Using the uniform exponential decay of L( 2 :^,ry^) ^ nj((z^, ry^) and (37), we finally end up with 


nj((^,ry)Wo = 0. 


Since W n belongs to the generalized eigenspace of L associated with those common eigenvalues of M(z, ry) 
and L(z, ry) in we can conclude that W n equals zero. Applying the induction relation (40), the whole 
sequence is zero which yields the expected contradiction. □ 

The crucial property that we use in the proof of Lemma 3 is the fact that up to z G , the eigenvalues 
of M(z,ry) lie either in D or . For the leap-frog scheme, this property would not be true if we had 
imposed the auxiliary numerical boundary condition + u” rather than 2 -|- A a ) • 

Let us also observe that we have used the fact that and a-r^ are nonconstant in order to study 
the induction relation (36). There might be some schemes for which Op^ and/or a-n are constant but for 
which one can still apply similar arguments as in the previous proof, even though (36) is an induction 
relation with fewer steps than (35). In this respect. Assumption 3 might be relaxed in specific applications. 
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Remark 2. The auxiliary problem (27) is in general not of the same form as (5) because in (27) one 
has to impose infinitely many numerical boundary conditions. This is due to the fact that the stencil 
of M incorporates points “on the left” with respect to the first space variable. A remarkable exception 
occurs for explicit schemes with s = 0, for in that case the multiplier M reads and (27) is exactly 
the auxiliary problem considered (and labeled (2.7)) in [CGll] where one imposes Dirichlet boundary 
conditions on finitely many boundary meshes (just use 5 ” = 0 for ji < —ri). The reader can then check 
that the energy-dissipation balance law of Proposition 2 in that case (s = 0, Qi = I) coincides exactly 
with the algebra involved in the derivation of the estimate (2.12) in [CGll]. The reader can also check 
that for s = 0, and Qi = I, Lemma 3 becomes a rather trivial exercise... 

There still remains the problem of constructing a set of dissipative numerical boundary conditions of 
the same form as (5) with s > 1, that is with finitely many numerical boundary conditions, and for which 
one can prove by hand both a semigroup and a trace estimate as in Theorem 2. 


3.3 End of the proof 

As explained in the introduction of Section 3, the linearity of (5) reduces the proof of Theorem 1 to the 
case {Fj') = 0, (< 7 "’) = 0, since we have already dealt with the case of zero initial data. We thus focus on 
(5) with (T”) = 0 and {g))) = 0, and write the corresponding solution (u”) as Fj = Vj" + w^, where the 
sequence (u”) solves: 

Lv'j = 0 , j € , ji >1, n>0, 

M u” = 0, j G Z'^, ji <0, n>0, (43) 

j G Z'^ , n = 0 ,..., s , 


= f 


and (wf) solves: 
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= 1 - ri,.. 


a=0 






0 

II 
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II 


(44) 


For Uj + w)[ to coincide with the solution (u”) to (5), it is sufficient to extend the initial data fj, ■ ■ ■, fj 
by zero for ji < —ri, which provides with the initial data in (43) on all Z'^, and to define the boundary 
source term in (44) by: 


s+i 


~n+s+l _„,n+s+l 

^1 ■ j 




(45) 


cr=0 


We can estimate the solution (vf) to (43) by applying Theorem 2. In particular, the trace estimate: 
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for Pi = max(pi, gi + 1) gives (recall the definition (45) of 
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We can apply Theorem 1 to the solution (w'^) to (44) because the initial data in (44) vanish. We get: 


sup e 2 III-u;"'III 2 _|- ^ 

n>0 .1 n,+oo 


-2771 At |||„,,n|||2 

^ III ^ nil— ri,+oo 


n>0 


Pi 




w 


,n II2 

il>- 


<C Y Ate-27riAt 




n>0 


jl=^-rl 


n>s+l 




< 


cE 

cr=0 


pa III 2 

' III 1—ri,+oo • 


Combining with the similar estimate provided by Theorem 2 for (u”), we end up with the expected 
estimate: 


sup +00 + — 

n>o . yAt + l 


llk"lll?-n,+oo 
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Pi 
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-)scx: 
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which completes the proof of Theorem 1. 


4 Conclusion and perspectives 

Let us first observe that in [Wad90], Wade has constructed symmetrizers for deriving stability esti¬ 
mates for multistep schemes, even in the case of variable coefficients. His conditions for constructing a 
symmetrizer are less restrictive than Assumption 2. However, the symmetrizer in [Wad90] is genuinely 
nonlocal and it is therefore not clear that it may be useful for boundary value problems. The main nov¬ 
elty here is to construct a local multiplier whose properties allow for the design of an auxiliary dissipative 
boundary value problem. This is our key to Theorem 1. 

In this article we have always discarded the dissipation term provided by the nonnegative form Dq . For 
the approximation of parabolic equations, this term may give some extra dissipation, but a crucial point 
to keep in mind is that the coefficients of the numerical scheme are assumed to be constant (which may 
in turn yield rather severe CFL conditions for implicit approximations of parabolic equations). Hence it 
does not seem very clear that our approach will yield stability estimates with “optimal” CFL conditions 
when approximating parabolic equations. This extension is left to further study in the future. 

The main possible improvement of Theorem 1 would consist of assuming that only the roots to (8) 
that lie on are simple. Here we have assumed that all the roots, including those in D are simple. If 
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we could manage to deal with multiple roots in D, then Theorem 1 would be applicable to numerical 
approximations of the transport equation (12) that are based on Adams-Bashforth methods of order 3 or 
higher (such methods have 0 as a root of multiplicity 2 or more at the zero frequency). 


The results in this paper achieve the proof of a ’’weak form” of the conjecture in [KW93] that strong 
stability, in the sense of Definition 1, implies semigroup stability. However, an even stronger assumption 
was made in [KW93], namely that the sole fulfillment of the interior estimate 
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when both the initial and boundary data for (5) vanish, does imply semigroup stability. The analogous 
conjecture for partial differential equations seems to be still open so far, but we do hope that our multiplier 
technique may yield some insight for dealing with the strong form of the conjecture in [KW93]. 
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